Implicit and explicit solvent models for the simulation of a single polymer chain in 
solution: Lattice Boltzmann vs Brownian dynamics 
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We present a comparative study of two computer simulation methods to obtain static and dynamic 
properties of dilute polymer solutions. The first approach is a recently established hybrid algorithm 
based upon dissipative coupling between Molecular Dynamics and lattice Boltzmann (LB), while the 
second is standard Brownian Dynamics (BD) with fluctuating hydrodynamic interactions. Applying 
these methods to the same physical system (a single polymer chain in a good solvent in thermal 
equilibrium) allows us to draw a detailed and quantitative comparison in terms of both accuracy 
and efficiency, ft is found that the static conformations of the LB model are distorted when the 
box length L is too small compared to the chain size. Furthermore, some dynamic properties of 
the LB model are subject to an L _1 finite size effect, while the BD model directly reproduces the 
asymptotic L — » oo behavior. Apart from these finite size effects, it is also found that in order to 
obtain the correct dynamic properties for the LB simulations, it is crucial to properly thermalize 
all the kinetic modes. Only in this case, the results are in excellent agreement with each other, as 
expected. Moreover, Brownian Dynamics is found to be much more efficient than lattice Boltzmann 
as long as the degree of polymerization is not excessively large. 



I. INTRODUCTION 

The rich variety of conformations which leads to many 
different intrinsic properties of polymer solutions has con- 
tinuously drawn considerable interest in soft matter re- 
search. Computer modeling is increasingly being used 
as an integral part of theoretical study, in order to both 
test existing theories and to trigger the development of 
new concepts. Furthermore, computer simulations have 
also become an essential tool in materials research, espe- 
cially for predicting and understanding the behavior of 
complex systems, where a complete theory is not avail- 
able. It has been proven to be an effective and inex- 
pensive way to study these systems. In order to observe 
large-scale properties, it is crucial to reduce the compu- 
tational cost by coarse-graining the details of the atomic 
structure. This is particularly true for polymer systems 
and studies of their universal static and dynamic prop- 
erties [H, 0. In this context, using a conventional bead- 
spring chain model to represent a polymer molecule in 
Molecular Dynamics (MD) simulations is usually suffi- 
cient 0, 0, H IE 13 ■ In the case of dilute and semidilute 
polymer solutions, a correct model also needs to take into 
account the effect of solvent molecules. This effect is two- 
fold: On the one hand, the good solvent quality results 
in swelling of the random coil; on the other, the solvent- 



mediated long-range dynamic correlations between dif- 
ferent segments of the chain, known as hydrodynamic 
interactions (HI), significantly influence the dynamical 
behavior [H, @, Q . In the present methodological study, 
we focus on the dilute regime which is theoretically most 
thoroughly understood. 

In order to capture hydrodynamic interactions in MD 
simulations, the solvent particles need to be incorporated 
explicitly. Typically, the number of solvent particles re- 
quired for such a model is of the order of thousands even 
for a short chain. Although such studies are feasible 
they are rather inefficient, for this reason. Therefore, a 
more coarse-grained description of the solvent is highly 
desirable. Two complementary approaches have been de- 
veloped to do this. "Mesoscopic" methods keep the sol- 
vent degrees of freedom, but describe them in a simpli- 
fied fashion. These include Dissipative Particle Dynamics 
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Multi-Particle Collision 
and lattice Boltzmann 
|2J,|25|, US H, Si. These 
approaches are typically one to two orders of magnitude 
faster than MD 1271. Conversel y, B rownian Dynamics 
(BD) simulations [3ClL l3ll l32l l33l l34j remove the solvent 
degrees of freedom completely, and take their effect into 
account via non-trivial long-range dynamic correlations 
in the stochastic displacements. This is possible due to 
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the time scale separation between the fast solvent motion 
and the slow conformational polymer degrees of freedom. 
Since the number of degrees of freedom is reduced drasti- 
cally, the method has the potential to save CPU time by 
additional several orders of magnitude, in particular in 
the dilute limit. However, a simple implementation of the 
correlations [3(| leads to an algorithm which scales like 
0(iV 3 ), where N is the number of Brownian particles, 
and therefore becomes infeasible as soon as N exceeds a 
few hundred [32| . It is therefore very important to treat 
HI by means of Fixman's algorithm [35[ (scaling roughly 
as 0(N 2 - 25 )), which we do in the present study. The 
recently introduced method by Geyer and Winter [36| 
would reduce the necessary CPU effort by roughly one 
additional order of magnitude for typical chain lengths, 
and exhibit a more favorable scaling (0(N 2 )). However, 
it is based upon an inexact approximation of the hydro- 
dynamic correlations that cannot be improved systemat- 
ically (in contrast to Fixman's method). For this reason, 
this algorithm was not implemented. 

None of these approaches is sufficient to reach N ~ 
10 3 . . . 10 4 ; this latter goal is only attainable by the im- 
plementation of very recent "superfast" BD alg orithms 
based upon Fast Fourier Transforms [13, [H, |39| . These 
latter algorithms scale as N 1+x log N, where the expo- 
nent x depends on the details of the underlying physics, 
and is usually substantially smaller than unity. These 
methods require the study of a confined system, and 
hence are not used in the present study. 

While the advantages and disadvantages of the meth- 
ods are well-known in general terms (and have resulted in 
differing methodological preferences in different groups of 
researchers), not much is known quantitatively in terms 
of a clear comparison of computational efficiency. The 
present paper aims at partly filling this gap. 

Recently, one of the present authors [2(| [27| has pro- 
posed a new mesoscopic method for simulating polymer- 
solvent systems. The solvent is represented by a fluid on a 
grid, simulated via the lattice Boltzmann approach, while 
the motion of the polymer chain is governed by a contin- 
uous MD model. The two parts are coupled by a simple 
dissipative force. The lattice Boltzmann (LB) method 
was originally developed to simulate hydrodynamics on 
a grid [Tjl, [22] • It has been shown that the LB method 
is a fast and effective method for simulating fluid flows, 
which has the sam e sp eed and accuracy as other Navier- 
Stokes solvers [H Hfl SJ S3 ■ Ladd [20, HI successfully 
applied the LB method to colloidal systems (originally 
with a conservative coupling) and showed that the CPU 
cost scales linearly with the number of particles. More- 
over, he showed how fluctuations can be incorporated 
into the LB model, which is essential in order to investi- 
gate Brownian motion [20 ]. This procedure has recently 
been refined and improved [H, [2!| . The dissipative cou- 
pling method [26|, [27| was thoroughly tested by applying 
it to a single polymer chain in solution, for which the 
data of a previous MD simulation [f| were available, and 
whose parameters were used as an input for the meso- 



scopic model. 

In this work, we study the dynamics of a single chain 
in a solvent to compare the predictions of the explicit sol- 
vent model via the LB method with the predictions of the 
implicit solvent model by BD simulations. Many recent 
simulation studies have investigated very similar systems, 
using various mesoscopic solvent models like MPCD [4lj . 
DPD [42j, direct solution of the Navier-Stokes equation 
[HI , and smoothed DPD [44[ . These studies were mainly 
done in order to test and verify the validity of the chosen 
mesoscopic simulation approach. Meanwhile it is fair to 
say that a single chain in solution is a standard bench- 
mark test bed for mesoscopic simulation methods, the 
necessary condition for passing being the correct repro- 
duction of the Zimm scaling laws within the limitations 
of finite chain length and finite solvent volume. The 
present study, however, aims directly at a quantititave 
comparison between two different methods, and there- 
fore it is crucial that that the underlying polymer model 
is exactly identical for both methods. This is the rea- 
son why the data of Refs. [HI E3, S H], though all 
exhibiting essentially the same physics, are of no direct 
relevance for the present investigation, since all of them 
employ slightly different polymer models, and simulate 
solvents with somewhat different viscosities. One very 
recent study by Ladd et al. [45[ has done a rather simi- 
lar comparison between LB and BD, and arrived at sim- 
ilar conclusions; the present paper should be viewed as a 
complementary study that puts some more emphasis on 
the issue of computational efficiency. 

After introducing the models, we show how to map 
the input parameters of the hybrid model onto the in- 
put values of the BD model to directly compare the pre- 
dicted quantities (Sec. [TTJ) . Section Hill then confirms the 
expected physical equivalence of the two approaches in 
terms of comparing static and dynamic data. Further- 
more, this section also presents our comparison on the 
numerical cost or benchmark data for both methods. Fi- 
nally, in Sec. IIVI we summarize the results and give our 
conclusions. 



II. MOLECULAR MODEL AND SIMULATION 
METHODS 

A. Molecular model 

In this work, a polymer molecule is represented by a 
conventional bead-spring chain model, which consists of 
N beads that are connected via N — 1 finitely extensible 
nonlinear elastic (FENE) massless springs. The Lennard- 
Jones potential, which acts between all monomers, is used 
to model the excluded volume (EV) effect. The two po- 
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tentials Vfene and Vlj are given by the expressions 

1/ ^FENE^o , A /" r \\ 

W = —^^-{wj J' 

^ = *(£-£ + i)' r<2V- ff , (2) 

where r is the bead-bead distance, /cfene is the spring 
constant and Rq the maximum extension of the bond, 
e and a are the energy and length parameters of the 
Lennard-Jones potential, respectively. 



B. The lattice Boltzmann (LB) method 

In this method, the evolution of the LB variables rij 
is g overned by the following lattice Boltzmann equation 



Hi (r + c 4 At, t + At) = m (r, t) (3) 

6 

+ E L n M r ' *) - n ?(P> u )) + <( r > *)• 
j'=i 

The variable ni(r,t) is the (partial) fluid mass density 
at grid site r at time t, corresponding to the discrete 
velocity Cj. At is the time step, and the lattice spac- 
ing is denoted by a. The small set of velocities 
(i = 1, . . . , 6, where the value of b depends on the de- 
tails of the model) is chosen such that c^Ar is a vector 
leading to the ith neighbor on the grid. Ljj is a collision 
operator for dissipation due to fluid particle collisions, 
such that the populations always relax toward the local 
pseudo-equilibrium distribution n^ q that depends on the 
local hydrodynamic variables p = ^ i u, (the total mass 
density) and u = J2i n i c i/ J2i n i ( tne local flow veloc- 
ity). The collision process is constructed in such a way 
that it conserves both p and u. n£(r, t) is the stochas- 
tic term, which is essential in order to simulate thermal 
fluctuations that drive Brownian motion. 

The local pseudo-equilibrium distribution can be rep- 
resented as a second-order expansion of the Maxwell- 
Boltzmann distribution, given by [23[ 
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n i q (p, u) = pw Ci 1 + 
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where w Ci are a set of weight factors, which depend 
on the sublattice i (i. e. the magnitude of Cj) and 
c s = \JTJZ{a/ At) is the speed of sound. In this work, 
we have used the algorithm proposed in Ref. [27| . how- 
ever, with the modification that the original 18-velocity 
model (D3Q18) was replaced by the D3Q19 19-velocity 
model [23| . The set of Cj consists of the particle be- 
ing at rest, the 6 nearest and 12 next-nearest neigh- 
bors on a simple cubic lattice. The magnitudes of the 
velocities corresponding to these three sets of particles 
are Cj = |cj| = 0, a/ At, and ^/2q/At, respectively. 



The weight factors for the D3Q19 model are wq = 1/3, 
toi = 1/18 and = 1/36. 

The early model of Ref. [27j only considered the ther- 
malization of modes related to the viscous stress tensor. 
It is important to note that even though this procedure 
is correct in the hydrodynamic limit, it provides poor 
thermalization on smaller length scales [28|. Adhikari et 
al. [28j have shown that by applying thermalization to 
all nonconserved modes one gets a significantly improved 
numerical behavior at short scales; the theoretical back- 
ground is now thoroughly understood [25l . |29(. In this 
work, we have also investigated the effects of thermaliza- 
tion of the kinetic modes on various dynamic properties. 

The coupling to the beads is done via simple interpola- 
tion of the flow velocity from the surrounding sites, and 
by introducing a phenomenological Stokes friction coeffi- 
cient (bare of the beads. This gives rise to a friction force 
on the particles, plus a Langevin force that balances the 
frictional losses. The total momentum is conserved by 
subtracting the corresponding momentum transfer from 
the surrounding fluid. It can be shown that this pro- 
cedure satisfies the fluctuation-dissipation theorem [25| |. 
For further technical details on this method and its the- 
oretical analysis, we refer the reader to Ref. (25j . 



C. Brownian dynamics (BD) simulations 

The configuration of a bead-spring chain is specified 
by the set of position vectors (i — 1,2, ... , N). The 
time evolution of this configuration is governed by the 
Ito stochastic differential equation [34|, 0, [i?} 



r,(t + At) = Ti(t) + (k B T) 1 Dy • (F; 
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J = 1,2,. 



(5) 



where summation over repeated indices is implied. Here 
the symbols At, fee, T, and denote the time step, 
Boltzmann's constant, the temperature, and the diffusion 
tensor, respectively, where the latter decribes the hydro- 
dynamic interactions between the beads. The forces F^- 
and F" 1 * are the spring and excluded-volume contribu- 
tions, repectively. W, are random variables representing 
a discretized Wiener process, such that (Wj) = and 



(W, g> W., 



ldij, where Sij 



is the Kronecker delta and 



1 is the unit tensor. Finally, the tensor By is related to 
the diffusion tensor such that Dy = • B J fc [4(| . 

The frictional properties of the chain and the hydro- 
dynamic interactions between the beads are modeled via 
the diffusion tensor D y . Its diagonal elements contain 
the bead friction coefficient ( = 6Tn] s d, where t] s is the 
solvent viscosity and d is the Stokes radius of the bead. 
The off-diagonal elements represent the hydrodynamic 
interactions via the tensor fly , for which we take the reg- 
ularized Rotne-Prager-Yamakawa (RPY) tensor [H, [4!| 
with solvent viscosity rj s and Stokes radius d. Taken to- 
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gether, the diffusion tensor is given by 



cMyi + (i-%) a 



(6) 



Further details can be found in Ref. [47[. 

The computationally most intensive part is to deter- 
mine the matrix By. Generally, Cholesky decomposi- 
tion of D.y is used to obtain By as an upper (or lower) 
triangular matrix, and the computational cost for this 
method scales as N 3 [35j]. Fixman made use of the fact 
that there are many possibilities to define By as some 
square-root matrix of Dy , and, noting that it is the vec- 
tor By • Wj that is required rather than the matrix By , 
applied a truncated Chebyshev polynomial expansion to 
obtain By • Wj with a lower computational cost, scaling 
roughly as iV 2 25 (35[. In the present paper, we make use 
of this accelerated technique as well. 



D. Unit systems and parameter mapping 

For implementation on a computer, physical quantities 
must be represented in certain units, i. e. in terms of 
suitable dimensionless ratios. This is typically done by 
choosing a natural unit system where three independent 
elementary quantities are set to unity. 

In the coupled MD/LB simulation approach, this is 
usually done by cho osing t he Lennard- Jones parameters 
e, cr, and t (t — ^Jmcr 2 /e, where m is the mass of the 
monomer) as the units of energy, length, and time, re- 
spectively; this choice enables one to make direct contact 
with MD simulations with explicit solvent |27| . Con- 
versely, BD simulations have traditionally [3J, [47J used 
a unit system where one chooses k^T as the energy unit 
and Ik = V^bTV^fene as the length unit. This is par- 
ticularly useful for the simulation of pure Gaussian (har- 
monic) chains where the interaction potential has neither 
an energy scale nor a length scale built in. The time unit 
in BD simulations is usually chosen as Tk = CA^fene)- 

Of course, a meaningful comparison of results requires 
that all data are represented in one common unit sys- 
tem. At this point, one realizes that this is less straight- 
forward than one might think at first glance. While the 
conversion of length and energy units is trivial, and di- 
rectly facilitated by the fact that both methods use the 
same molecular model for the polymer chain, and per- 
form the simulations at the same temperature, the con- 
version of time units is not. This is so because of the 
different time scales underlying the basic updating algo- 
rithms: The MD/LB method is based upon simulating 
the system on inertial time scales, while BD focuses di- 
rectly on the larger diffusive time scales. This is directly 
reflected by the occurence of the intertial parameter m 
(monomer mass) in the MD/LB time unit, which does 
occur in the BD model, and the diffusive parameter C 
(monomer friction constant) in the BD time unit. It is 
important to notice that C is not an input parameter to 



the MD/LB model. Rather, one needs to carefully distin- 
guish between the short-time friction coefficient Chare that 
is indeed an input parameter — it describes the Stokes 
coupling of the monomer to the LB fluid in its immedi- 
ate vicinity — and the long-time friction coefficient C e ff 
describing the particle's long-time response that is modi- 
fied by solvent backflow effects. It is this latter parameter 
that must be identified with the BD £, and it is esentially 
an output parameter. Fortunately, the relation between 
Char e and Ceff is well understood — Ahlrichs and Diinweg 
I25L l27l have shown that 
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(7) 



where g is a numerical prefactor and a is some measure 
for the range of interpolation to the surrounding lattice 
sites. For linear interpolation to the nearest sites, one 
finds g w 25, if a is the LB lattice spacing. For highly ac- 
curate results, one should also take into account a small 
correction for the finite size of the simulation box [25| : 
this has however not been done in the present pa per . 
Rather, we took the mapping determined in Ref. [271 ] 
to calculate Ceff from Chare and identified this with the C 
parameter of the BD calculations. Although the phys- 
ical mapping done in this way is essentially correct, it 
is important to notice that this aspect introduces a cer- 
tain amount of numerical inaccuracy when it comes to 
quantitative comparisons. 

Given the fact that it is intrinsically impossible to run 
the two simulations with the same unit system, we chose 
to keep the previously established systems of the two re- 
spective methods, and to map the results a posteriori 
by the procedure outlined above. Furthermore, we chose 
to present all results in MD/LB units. Technically, this 
means that for length and time unit conversions we need 
to set la = l*lk and It = t*Tk, where the "*" superscript 
denotes BD non-dimensionalization, while "-" denotes a 
non-dimcnsionalization for MD/LB. The corresponding 
factors for the conversion from BD to MD/LB units are 
then trivially found to be 
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(8) 
(9) 



note that /cfene and Ceff are &fene and Ceff in non- 
dimensional MD/LB units, respectively, while the ra- 
tio ksT/e is just the non-dimensionalized temperature 
in MD/LB units. 

The dimensionless hydrodynamic interaction parame- 
ter h* used in the BD simulations is essentially a non- 
dimcnsionalized Stokes radius. We thus find 



h* ' \fiilk = d 



Ceff CcffC 



(10) 
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or 

We therefore parameterized our simulations by first 
picking simulation parameters for the MD/LB model 
(which were then directly used for MD/LB), then con- 
verting these to BD units using the procedure outlined 
above, and then running the thus-obtained equivalent 
BD model. For these latter simulations, a time step size 
At* = 0.005 (in BD units) was found to produce suffi- 
ciently accurate results. 

E. Choice of parameters 

The physical input values for the present model are 
chosen from the benchmark values developed in Ref. [27[ , 
which have been shown to reproduce the results of a typ- 
ical pure MD simulation (fj. As in the comparison be- 
tween LB and MD simulations, we study a system of 
a single polymer chain of length N = 32 monomers im- 
mersed in a fluid with temperature k^T/e — 1.2, density 
p = 0.864, and kinematic viscosity v = 2.8. The lattice 
spacing a is set to unity, which is roughly identical to the 
bond length; this is necessary to resolve the hydrody- 
namic interactions on small length scales with sufficient 
accuracy. 

Furthermore, following Ref. [27j], the coupling param- 
eter Cbare was set to 20.8. The values of the FENE spring 
potential parameters are /cfene = 7 and Rq = 2. The 
time step size for the polymer (the MD part of the sim- 
ulations) is set to At = 0.01. It should be noted that 
such a large time step is possible since the inclusion of 
dissipation and noise lead to a substantial stabilization, 
compared to purely microcanonical MD. The value of 
the time step that updates the fluid should be chosen 
in a way such that the LB variables rij do not become 
negative too often. Here, we choose Af = 0.02, where 
such a case rarely occurred during the observation time. 
It is important to mention another free input parame- 
ter which governs the time scale for the evolution of hy- 
drodynamic interactions, known as the Schmidt number 
Sc = is /Do, where Do is the diffusion constant of the 
single monomer. This parameter can be set arbitrarily 
in the LB method by choosing is and Dq (which can be 
tuned by choosing (bare) accordingly. Ideally, the value 
of Sc should be chosen such that hydrodynamic interac- 
tion evolves much faster than the diffusion of a monomer. 
In our case, we have Sc ss 32, which has been shown to 
result in Zimm-like behavior [53, |42j ■ 

For the LB simulations, the polymer chain moves 
within a cubic box of length L with periodic boundary 
conditions, while it is drifting freely in an infinite medium 
for the BD simulations. In order to accurately compare 
various properties between the two systems, one must 
understand the effects of the box length L on any ob- 
servable of interest in the LB simulations. Thus it is 



essential to only compare quantities under identical con- 
ditions (i. e. independent of the box length). Hence var- 
ious box lengths L ranging from 10 to 35 Lennard- Jones 
units were investigated; this allows us to extrapolate to 
the L — > oo limit. 



III. RESULTS AND DISCUSSION 

A. Static properties 

The mean square radius of gyration and the mean 
square end-to-end distance are given by 

( R D = a** £<»■&>> ( 12 ) 
(Rl) = ((rN-nf) (13) 

with r.jj = |r.j — rj\ being the inter-particle distance. 

These two quantities are both related to the number 
of monomers by the expression 

(Rl) ex (Rl) ex N 2u , (14) 

where is is the Flory exponent. For a self-avoiding walk 
(SAW), the Flory exponent v is 0.588 [50j. In principle, is 
can be obtained from simulations, using the scaling law in 
Eq. 1141 However, this method would require simulations 
for a wide range of N values. Alternatively, one can use 
the static structure factor 

S W = -^X;<«*p(ik.r«)) 

^3 

to obtain v much more efficiently 

In the scaling regime R~ x <C k <C ao (do being a mi- 
croscopic length of the order of the bond length), a power 
law relation between the static structure factor and the 
wave vector k holds: 

S(k) oc k" 1 ^. (16) 

Figure [T] shows the static structure factor as a function 
of wave vector k for the LB simulations with the presence 
of thermalization of all modes, and the BD simulations. 
It can be clearly seen that the values of the static struc- 
ture factor obtained from the LB simulations are exactly 
the same as those obtained from the BD simulations, in- 
dicating that they have the same static conformations. 
From Eq. [T|)l the value of is can be extracted from the 
linear region of the log-log plot of S(k) vs k. As ex- 
pected, the values for is obtained via this method are the 
same for both the LB and the BD simulations, as re- 
ported in Table [U However, they are approximately 5% 
higher than the asymptotically correct value, which is a 
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consequence of the finite chain length. The results for 
the mean square radius of gyration and the mean square 
end-to-end distance in Table Q] further confirm this agree- 
ment with regard to static conformations between the 
two methods. However, at small box length [L = 10), 
the results for these static properties for the LB method 
deviate from their asymptotic values. The discrepancy 
observed here always arises when the box length is too 
small compared to the chain size, where the chain is more 
likely to wrap over itself due to spatial restriction and 
hence alter its static conformations. We also found that 
the two versions of LB thermalization ("stresses-only" 
vs "full" thermalization) yield identical results for the 
chain conformational statistics. In general, we only quote 
values obtained for full thermalization, unless indicated 
otherwise. 

The hydrodynamic radius for a single chain in an infi- 
nite medium is given by 




For a chain in a finite box, as is the case here in the LB 
method, it has been shown that the hydrodynamic inter- 
actions of the chain with its periodic images effectively 
increases i?n H, H3 • In order to account for this finite- 
size effect, a finite-size correction of order L~ 1 for most 
dynamic properties, resulting from the slow r _1 decay of 
hydrodynamic interactions, is required @, The re- 
sults for the infinite-box value (-Rjj ) agree excellently 
with each other for all simulations (see Table [I]) . Since 
the overwrapping effect is more sensitive to large inter- 
particle distances, it turns out that the deviation in the 
inverse hydrodynamic radius is too small, for the range 
of box lengths used, for it to be distinguishable. As can 
be seen in Table HI the deviation is more pronounced for 
the radius of gyration, and even more for the end-to-end 
distance. 



B. Dynamic properties 

According to dynamic scaling, the longest relaxation 
time tz of the chain is, by order of magnitude, identical 
to the time that the chain needs to move its own size, i. e. 
DcmTz ~ Rg, where Dcu is the diffusion constant of the 
chain's center of mass. This leads to a dynamic scaling 
law rz cx iig, where z is the dynamic scaling exponent. 
For a chain with hydrodynamic interactions, this relax- 
ation time is known as the Zimm time rz. For this case, 
Dcm cx R^ 1 in the limit of long chains. This implies that 
rz cx i?g, which gives a dynamic exponent of z — 3 for 
models with HI. For the Rouse model (i. e. chains with- 
out hydrodynamic interactions), where Dqm cx N , one 
finds a dynamic exponent of z = 2 + \ jv. These quanti- 
ties will be referred to in the discussion below. 

The mean-square displacement of the chain's center of 



mass 

g 3 (t) = ((Rcm(*o + 1) - RcM(io)) 2 ) (18) 

for both methods is depicted in Fig. [51 From the fig- 
ure, it can be clearly seen that 53 strongly depends on 
the box length L for the LB simulations. Moreover, they 
seem to converge to the value predicted by the BD sim- 
ulations [L — 00) in the limit of large L. Effects of 
thermalization of the kinetic modes in LB simulations on 
this property will be discussed subsequently. The chain's 
center of mass diffusion constant Dqm can be determined 
by the slope of the gs vs t curve, where the relationship 
gs{t) — 6DcMt holds. By fitting a power law to the sim- 
ulation data, we obtain the exponents and the diffusion 
constants shown in Table |TJ These exponents support 
the prediction of simple diffusive behavior (i 1 ). Theoret- 
ically, one would expect that two diffusive regimes exist: 
On the one hand, there should be a short-time diffusive 
regime, corresponding to time scales well below the Zimm 
time, t < tz, but also well above the ballistic regime, 
t ^ To; note that To > only in the LB case, since the 
BD equation of motion is overdamped. On the other 
hand, there should be free diffusion for times t > rz. 
Both these regimes exhibit t 1 behavior, but with differ- 
ent prefactors, with a smooth crossover around the Zimm 
time [HI, HH, HH • In principle these two different diffusion 
constants can be obtained via fits to the corresponding 
regimes. In practice, however, it turns out that the val- 
ues are very close to each other, and hence the crossover 
is very smooth [32|, [5l], . Therefore its unambiguous 
identification is very difficult, i. e. impossible within the 
resolution of our data. 

The mean-square displacement of a single monomer i 
is given by 

gi(t) = ((ri(t +t)- r .-(*o)) 2 ). (19) 

Here, only the two innermost monomers near the center 
of the chain are evaluated to eliminate end effects; the 
results are plotted in Fig. [3J The values of gi behave 
similarly to those of g 3 . In the sub-diffusive time regime, 
corresponding to the short-time diffusive regime for 173, 
here evaluated between i = 20 and 80, the scaling behav- 
ior gi (t) cx t 2 l z is predicted [|J ■ The corresponding expo- 
nents obtained from a power-law fit are listed in Table U 
and indicate a value of z — 2.75 as L — > 00. Regardless 
of the finite-size corrections due to the box length and 
the effects of thermalization, these values clearly favor 
the Zimm model compared to the Rouse model, which 
predicts gi(t) cx £ a54 . The deviation from the asymp- 
totic Zimm value z = 3 (or g\ cx i 2 / 3 ) is mainly a result 
of finite chain length. 

Figure 3] shows the mean-square displacement of a sin- 
gle monomer in the center of mass system (i. c. the two 
innermost monomers to eliminate end effects) 

g 2 (t) = (([r < (io+t)-R OM (io+*)] (20) 
- [ri(i )-RcM(io)]) 2 >- 
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Interestingly, when viewed within the center of mass 
system, all the results lie on top of each other, regard- 
less of the box length L. This result also holds for LB 
simulations without full thermalization. This shows that 
the global center-of-mass motion of the chain is actually 
the primary contribution to the deviations between LB 
and BD results. In the data of Fig. [4] this contribution 
is suppressed: In terms of Rouse modes, only the in- 
ternal modes remain. For these modes, however, it has 
been shown [23] that the HI with the periodic images 
is much weaker, while the leading-order HI cancels 
out. Therefore, the corresponding finite-size effect scales 
as L~ 3 instead of and this is so small that it is 

invisible in Fig. |U 

Theoretically, these data can also be used for estimat- 
ing the Zimm time as the time where the crossover to the 
long-time plateau occurs. However, the crossover is quite 
extended and smooth, making it difficult to extract. We 
therefore estimated the Zimm time via 



Ttr 



6L>cm 



(21) 



Strictly speaking, this definition is only valid for a sin- 
gle chain in an infinite medium, where there is no finite 
box size effect. In the presence of finite box size, it be- 
comes the definition for the translational time (r tr ) rather 
than the Zimm time: The former is subject to an L^ 1 
finite-size effect, due to the strong L-dependence of -Dcm, 
while the latter, being defined via the relaxation of in- 
ternal modes, is only subject to an L~ 3 size effect, as 
discussed above. The translational times obtained from 
Eq. [21] (as shown in Table fl} are indeed different for dif- 
ferent box lengths L, as expected. Conversely, the results 
displayed in Fig. |4] indicate that the systems with differ- 
ent box sizes have (essentially) all the same (internal- 
mode) Zimm time, since their data all lie on top of each 
other. 

Next, we focus on the leading order L^ 1 finite size cor- 
rection for the long-time diffusion constant of the chain's 
center of mass, -Dcm- In principle, a plot of Dqm vs 
should give a straight line for large L, and an extrapo- 
lation to the limit L — > oo should yield the same value 
as predicted by the BD simulations. Figure [5] shows the 
values of Dqm for the LB simulations with and with- 
out thermalization of the kinetic modes at various box 
lengths L plotted together with the value obtained from 
the BD simulation at L = oo. It is worth mentioning 
that the BD value of Dqm can be obtained from the 
mean square displacement of the chain center of mass 
or via Fixman's expression (5l| . The latter method has 
been shown to produce a much more reliable result and 
is easier to carry out [32l ]. The value reported here has 
been cross checked by both methods and the results are 
almost the same within error bars. For the LB simula- 
tions without thermalization of the kinetic modes, the 
value of Dcm at the asymptotic limit L = 00 is differ- 
ent from that predicted by the BD simulations by about 
9.5%. However, when all the kinetic modes in the LB 



simulations have been thermalized, the deviation in Dcm 
reduces to 3%. This result clearly indicates that it is very 
important to thermalize all the kinetic modes in order to 
obtain correct values for dynamic properties. 

The reason for the remaining small discrepancy be- 
tween LB and BD is not completely clear, since there are 
numerous possible sources. Firstly, it should be noted 
that the underlying equations of motion are quite dif- 
ferent: LB works with inertia, while BD employs over- 
damped dynamics. This results in different Schmidt 
numbers Sc and different Mach numbers Ma, the latter 
being defined as the ratio of the flow velocity to the speed 
of sound: Both are finite in the LB method, while in the 
BD case they are strictly infinite (Sc) and zero (Ma), 
respectively. Furthermore, the shape of the HI function 
at small interparticle distances is somewhat different for 
the two methods: In the BD case, we employ the RPY 
tensor, while the nearest-neighbor interpolation for LB 
results in a short-range HI that differs somewhat from 
the RPY tensor (see also the discussion in Ref. HEj])- Fi- 
nally, it should be noted that the value of the constant g 
in Eq. which is crucial for the mapping between the LB 
friction parameter (bare an d the BD friction £ e ff, is only 
known with some numerical inaccuracy. For highly accu- 
rate mappings, it is also necessary to include a finite-size 
correction in the definition of g [25j ; this was not done in 
the present study. 

In order to examine whether the thermalization of the 
LB kinetic modes is also important for the internal modes 
of chain motion, we have performed a Rouse mode anal- 
ysis. The Rouse modes for a discrete chain are defined 
as 27, 53] 



1 



~N 



(22) 



for p = 1,2,..., N - 1. 

Within the approximation of the Zimm model, the au- 
tocorrelation function of the modes should decay expo- 
nentially 



(X p (to + t) ■Xp(t)) 



= <'M> ( - — 

Tp 



(23) 



where r p is the relaxation of the p-th mode. To vali- 
date our Rouse mode analysis routine, we have carried 
out extensive simulations for a (Gaussian) Rouse chain 
of N — 8 in the absence of HI and EV; the results for 
t p are in excellent agreement with the analytical predic- 
tions [2| . Figure [5] shows the normalized autocorrelation 
function for p — 1, 2, . . . , 5 for the LB model with box 
length L = 25, and the BD model. For nonzero times, 
there is a small deviation between LB and BD, the latter 
exhibiting again a slightly faster dynamics. This devia- 
tion systematically becomes smaller upon increasing the 
mode index p. Since high mode index means essentially 
relaxation on a rather small length scale, it is tempting 
to attribute the deviation to the finite propagation of HI 
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in the LB model, i. e. to retardation effects, or effects 
of finite Schmidt number, which are more important on 
large length scales than on small ones. Nevertheless, this 
hypothesis is not proven. 

In Ref. [27| it was shown that the autocorrelation 
function is only subject to an L~ 3 finite size effect, in 
contrast to the usual L~ x behavior. Figure [7] shows the 
value of the autocorrelation function of the first Rouse 
mode Xi at a fixed finite time i — 700, for LB simu- 
lations at various box lengths L and BD simulations at 
L = oo. Within our numerical resolution, the data in- 
deed confirm this L~ 3 finite size effect, both with and 
without thermalization of the kinetic modes. Further- 
more, they demonstrate again that thermalization of all 
the kinetic modes in LB simulations improves the accu- 
racy of the dynamic properties and brings them closer 
to the BD prediction: The deviation in the extrapolated 
limit L — ► oo is reduced from 3% down to 2%. The rea- 
sons for the remaining discrepancies are probably of the 
same nature as in the case of Dcu ■ 

We have also evaluated the dynamic structure factor, 
which is defined as 

S ^ *) = E < ex P ( ik ■ - T M))) ■ ( 24 ) 

When both wave number and time are in the scaling 
regime (i. e. R~ 1 < k < 1 and t « ( « t z ), S(k,t) 
is predicted Q to exhibit the scaling behavior 

S(k,t) = S(k,0)f(k z t). (25) 

A plot of S{k,t)k,y y against (k z t) 2 l z should collapse 
to a single curve [23]. The results for both methods are 
shown in Fig. [5] The data were restricted to the scal- 
ing regime 20 < t < 80 and 0.7 < k < 1.5. These ranges 
were obtained from the single monomer mean-square dis- 
placement, Fig. [H and from the static structure factor, 
Fig. [TJ respectively. Here, we particularly focus on ad- 
justing the exponent z such that it would produce the 
best total data collapse for a chain in an infinite medium 
(i. e. in the BD model). Obviously, the results from the 
simulations show Zimm-like rather than Rouse-like be- 
havior. Even though we have suppressed the finite box 
size effect, a dynamic exponent of z = 2.75 yields the best 
data collapse, which is somewhat smaller than the cor- 
rect asymptotic one. This result is also consistent with 
the value of z obtained earlier via the exponent of g\ in 
the sub-diffusive scaling regime (i. e. 2/z — 0.728). The 
deviation from the asymptotic value is due to the finite 
chain size used here, and one can expect z — 3 only in 
the long chain limit N — > oo. 

More detailed comparisons of the structure factor 
S(k, t) are shown in Fig. [5] (k dependence at constant 
time) and Fig. [TU] (time dependence for the normalized 
structure factor S(k, t)/S(k, 0) at constant k). 

Figure [9] shows the structure factor for BD simula- 
tions for a wide range of k at three different times and 
the data clearly indicate that the structure factor de- 
cays rapidly with time. The normalized structure factor 



S(k,t)/ S(k,0) for three different k values are shown in 
Fig. [TU1 and the data seem to indicate that the LB results 
approach the BD data as L is increased, as expected. 

C. Efficiency 

For the ultra-dilute system considered here, the lattice 
Boltzmann part of the hybrid LB method uses up most 
of the computational resources as the CPU cost for the 
MD part for the polymer chain is negligible. Since the 
dynamic properties predicted by the LB model arc sub- 
ject to a finite-size correction of order L _1 , extrapolation 
is required to obtain these properties in the asymptotic 
limit L — ► oo. To perform this extrapolation, together 
with checking that indeed the asymptotic behavior 
has been reached, one needs the results of at least three 
different box lengths. Moreover, the box length should 
be large enough compared to the chain size such that it 
does not alter the static properties. The data displayed 
in Table [T] indicate that it is safe to choose L such that 
y/{B%)/L < 0.5. The three different box lengths L cho- 
sen here are yj (R%)/L = 0.5,0.4, and 0.3. In this work, 
we set the total CPU time required for the LB simula- 
tions to be the sum of all the CPU times required to run 
1000 MD time steps for each of the chosen box lengths. 
For BD, we take the CPU time needed to observe the 
system for the same time span in physical units. Each 
of the simulations performed for the CPU time compari- 
son was run on an Itanium 2 processor of a 1.6 GHz SGI 
Altix server 3700. All the parameters used to carry out 
this comparison are the optimal values for both meth- 
ods. Several chain sizes ranging from N = 16 to 1024 
have been used to obtain the CPU cost for comparison. 
The results are shown in Fig. QT] For the LB method, it 
is clear that the CPU cost scales linearly with the num- 
ber of particles, i. e. the number of grid point s that the 
solvent lives on, or L 3 . Since the ratio */ (R 2 ) / L is kept 
constant, or L oc yj (R 2 ) oc N u , this leads to a CPU cost 
scaling as N 3v . This is indeed found in our benchmarks, 
see Fig. 1111 Similarly, our data also confirm the pre- 
dicted N 2 - 25 CPU cost scaling for BD. Though the LB 
exponent is lower than BD, the large prefactor ensures 
that the total CPU cost for LB is much more expensive 
compared to BD for the typical chain lengths used in the 
literature. It is only when the chain length is excessively 
large (i. e. N of the order of 10 6 or higher) that LB will 
become superior to BD for a single-chain system. 

The situation completely changes if one studies a semi- 
dilute solution instead, as has been done in Ref. [HiJ . For 
such a system, we have not done a comparison between 
LB and BD in terms of actual simulations; however, by 
means of scaling considerations one can roughly estimate 
what the likely result of such a comparison would be. A 
scmidilutc solution comprises M chains of N monomers 
each, such that the total number of monomers is MN. 
Therefore the BD CPU cost scales as (MN) 2 - 25 , while 
the LB CPU cost depends on the density. Within the 
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blob picture of semidilute solutions, one views a chain 
as a sequence of "blobs" , each comprising n monomers, 
and having size £, which can be viewed as the typical 
correlation length of density fluctuations, or the typ- 
ical distance from which point on chain-chain interac- 
tions become important. Since the conformation statis- 
tics within the blob is that of a self-avoiding walk, one 
has £ ~ an v , where a is the monomer size. The sequence 
of blobs forms a random walk, hence R e ~ S,(N /n) 1 / 2 . 
This gives the minimum size of the simulation box, i. e. 
L ~ ^(N/n) 1 ' 2 ~ an v (N/n) 1 l 2 = aN^^/Nf- 1 ' 2 , or 
L 3 - a 3 N 3u {n/N) 3u - 3 / 2 . We thus see that the CPU ef- 
fort for the LB method is even slightly decreased by the 
factor (n/N) 3v ~ 3 ' 2 compared to the single-chain case at 
the same N, due to the shrinkage of the chains result- 
ing from excluded-volume screening. In order to esti- 
mate the number of chains M, we note that the arrange- 
ment of blobs is space-filling, i. e. L 3 ~ £ t 3 M(N/n) ~ 
a 3 Mn 3v (N/n) . Comparing this with the previous expres- 
sion for L 3 , one finds M ~ {N/n) 1 / 2 . Therefore the BD 
effort, compared to the single-chain case, is increased by a 
factor of M 2 - 25 ~ (N/n) 1 - 125 . Taken together, this means 
that the ratio between LB effort and BD effort is changed 
by a factor of - (jv/n) 3 " +1 - 125-1 - 5 ~ (N/n) 1A25 in favor 
of LB. For N/n — 30, which is needed as a minimum to 
resolve the Gaussian statistics of the chains as a whole, 
one obtains a factor of 130, which more or less compen- 
sates the two orders of magnitude seen in Fig. [11] Taking 
into account that for such a system the BD simulation 
would have to calculate the HI with the periodic images, 
e. g., via Ewald sums, which is much more complicated 
than the present single-chain simulation, one sees that 
for a semidilute solution clearly LB is more efficient, un- 
less a "superfast" BD algorithm [3?], [H, Hi| is used. For 
the latter case, the answer is not yet known. The results 
of Ref. [HI indicate that LB/MD may be favorable for 
a rather small number of monomers; however, this study 
was done (i) in a non-trivial geometry, which implies a 
more complicated BD method, and (ii) under complete 
neglect of thermal fluctuations in the LB simulations, re- 
sulting in a substantial reduction in CPU effort. Such 
a (partial or complete) neglect of thermal noise is some- 
times justified in strong nonequilibrium situations such 
as studied in Ref. [55j |; in that particular case, the jus- 
tification was checked by additional tests [56]. Another 
possible situation where LB noise is negligible is the case 
of strong coarse-graining, where a single lattice site can 
already be considered as a macroscopic thermodynamic 
system (for a detailed discussion, see Ref. ; 25]). However, 



in the general case, and certainly in thermal equilibrium 
or weak nonequilibrium, the proper inclusion of thermal 
noise is necessary, as demonstrated theoretically in detail 
in Ref. [25[ , and also corroborated by the present numer- 
ical results. For the general case, the estimate of the LB 
CPU effort given in Ref. [55[ is therefore too optimistic. 
IV. CONCLUSIONS 

The present study has shown that Brownian dynamics 
simulations are capable of reproducing various properties 
predicted by a hybrid LB/MD model (or vice versa). We 
have demonstrated how to obtain the input values for the 
BD simulations from the physical input parameters of the 
LB model such that both models would produce the same 
static and dynamic properties. For the LB model, most 
dynamic properties are subject to a finite-size correction 
of order L . In addition to this, it is very important 
to thermalize all the kinetic modes in order to obtain 
the correct dynamic properties. Those results that are 
not affected by L^ 1 finite size effects, such as the mean 
square displacement in the center of mass system, or the 
Rouse mode autocorrelation function, agree very favor- 
ably with each other. For highly dilute systems where the 
simulation of a single chain is sufficient, BD is usually the 
method of choice, as it is much more efficient than the 
coupled LB/MD approach, and finite box size effects are 
absent. The situation changes however in the semidilute 
case, where it is easy to estimate that BD will not be 
able to compete, unless "superfast" algorithms are used. 
Moreover, one should take into account that the hybrid 
LB/MD algorithm is rather easily adaptable to compli- 
cated boundary conditions, and can even be applied to 
flows at high Reynolds numbers, where the fluid degrees 
of freedom become intrinsically important, and cannot 
be handled in terms of a Green's function. 
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Tables 





LB 


BD 


Box length L 


10 


15 


25 


oo 


Time step 


0.02 


0.02 


0.02 


0.005 


exponent v 


0.615 ± 0.005 


0.617 ± 0.005 


0.619 ± 0.005 


0.619 ± 0.004 


(Re) 

(Rl) 


94.56 ± 1.20 


100.05 ± 1.26 


100.20 ± 1.28 


99.22 ± 1.24 


14.83 ± 0.10 


15.31 ± 0.11 


15.36 ± 0.11 


15.25 ± 0.11 




0.291 ± 0.0005 


0.290 ± 0.0005 


0.289 ± 0.0005 


0.290 ± 0.0005 


pi-exp. a 


0.640 ± 0.0005 


0.675 ± 0.0005 


0.710 ± 0.0005 


0.728 ± 0.0006 


pi-exp. a ' b 


0.645 ± 0.0006 


0.684 ± 0.0006 


0.714 ± 0.0006 


0.728 ± 0.0006 


<73-exp. a 


1.008 ± 0.0008 


1.020 ± 0.0008 


1.050 ± 0.0008 


0.995 ± 0.0008 


Dcm x 10~ 3 


3.914 


5.162 


6.959 


9.843 




±1 x 10~ 3 


±1 x 10~ 3 


±2 x 10~ 3 


±1 x 10~ 2 


ftr (estimate) 


631.36 ± 4.43 


492.01 ± 3.56 


368.51 ± 2.67 


258.28 ± 1.87 



TABLE I: Properties for a single chain of length N = 32 obtained from lattice Boltzmann simulations at various finite box 
lengths and Brownian dynamics simulations in infinite medium. a Exponent obtained by fitting a power law in the sub-diffusive 
scaling regime of the chain in Lattice Boltzmann simulations, i £ [20 : 80]. b Exponent obtained from Lattice Boltzmann 
simulations without thermalization of all the kinetic modes. 
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FIG. 1: The static structure factor for the LB simulations (at various box lengths L) and the BD simulations (L = oo) for a 
wide range of dimensionless wave vectors k. 




FIG. 3: The dimensionless mean-square displacement gi(t) of the central monomer, Eq.[l9] Values of the exponent z at various 
box length L in the sub-diffusive scaling regime are also listed in Table [I] 
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FIG. 5: The dimensionless long time diffusion constant for the center of mass at various box lengths L. 
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FIG. 6: Normalized autocorrelation function of the first 5 Rouse modes X p (Eq. I23[l for LB simulations at L = 25 and BD 
simulations at L —* oo. 




FIG. 7: The autocorrelation function for the first Rouse mode Xi at a finite time value of t = 700 for LB simulations at various 
box lengths L and BD simulations at L — > oo. 
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FIG. 8: Scaling plot of the dynamic structure factor for (a) Rouse scaling (z — 3.7, top), (b) asymptotic Zimm scaling (z = 3, 
center), and (c) z — 2.75 (bottom), which produces the best collapse. 
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FIG. 10: Time evolution of the normalized dynamic structure factor at three different k values. Data for k = 0.047T are displayed 
in the inset. 
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System size (N) 



FIG. 11: Comparison of the CPU time required by the LB and BD systems for the equivalent of 1000 time steps for a wide 
range of system sizes (chain lengths) N. 



